Carga de librerías y datos

library(ggplot2)
library(plotly)
library(tseries)
library(MASS)
library(forecast)
library(lmtest)
library(caschrono)
library(lubridate)
library(timeDate)
library(tsoutliers)
library(dplyr)

# Se limpia el espacio de trabajo:

rm(list = ls())

# Se asigna el directorio donde están los datos

# setwd("C:/Users/Javier/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")
setwd("C:/Users/jherraez/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")

datos <- read.csv("_11_CACERES.csv")

Dividimos en train y en set y se convierten los datos en un objeto de tipo time series para poder aplicar sobre ellos las funciones de R.

datos$FECHA = as.Date(datos$FECHA,format='%d/%m/%Y')
datos.train <- subset(datos, FECHA<=as.Date('01/12/2018',format='%d/%m/%Y'))
datos.test <- subset(datos, FECHA>as.Date('01/12/2018',format='%d/%m/%Y'))

datos.train.ts <- as.ts(datos.train$TARGET, frequency = 12)
datos.test.ts <- as.ts(datos.test$TARGET, frequency=12)

Representación gráfica de los datos

names(datos.train)
## [1] "FECHA"  "TARGET"
graficoInicial <- ggplot(aes(x= FECHA, y = TARGET), data = datos.train) +
  geom_line(color = '#d84519', size = 1) + 
  xlab('FECHA') + ylab('Matriculaciones')

ggplotly(graficoInicial)

Estacionalidad

Estacionalidad en varianza

Evaluamos si es necesario transformar la serie para hacerla estacionaria en varianza.

box_cox <- boxcox(TARGET ~ FECHA,
                  data = datos.train,
                  lambda = c(0, 0.5, 1))

lambda <- box_cox$x[which.max(box_cox$y)]
lambda
## [1] 0.1616162

Como lambda toma un valor cercano a cero, tomamos logaritmos.

datos.train$log_target=log(datos.train$TARGET)
datos.test$log_target=log(datos.test$TARGET)
datos.train.ts <- as.ts(datos.train$log_target, frequency = 12)
datos.test.ts <- as.ts(datos.test$log_target, frequency=12)

Estacionalidad en media

Probamos a hacer el test de Dickey-Fuller aunque no confiamos demasiado en su resultado, ya que es poco exigente.

adf.test(datos.train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  datos.train.ts
## Dickey-Fuller = -1.6659, Lag order = 5, p-value = 0.7159
## alternative hypothesis: stationary

Todo parece indicar que va a ser necesaria una diferencia regular, pero nosotros ajustamos directamente sin nada.

Ajuste de un modelo SARIMA a la serie

Observamos los gráficos f.a.s. y la f.a.p. para realizar el primer ajuste.

acf(datos.train.ts, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(datos.train.ts, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

Lo habitual es empezar con una estructura AR(1). La fuerte correlación de orden 1 en la f.a.p parece reforzar esta decisión, por lo tanto, se propone un SARIMA(1,0,0)x(0,0,0)12 + MU

ajuste1 <- Arima(datos.train.ts,
                 order = c(1,0,0),
                 seasonal = list(order = c(0,0,0), period = 12),
                 method = "ML")
ajuste1
## Series: datos.train.ts 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.9059  5.5001
## s.e.  0.0303  0.2199
## 
## sigma^2 = 0.09146:  log likelihood = -42.67
## AIC=91.34   AICc=91.47   BIC=101.11
coeftest(ajuste1)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       0.905855   0.030348  29.849 < 2.2e-16 ***
## intercept 5.500117   0.219941  25.007 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajuste1) # debajo de 0.8 todo, no correlaciones
##                   ar1   intercept
## ar1        1.00000000 -0.02691084
## intercept -0.02691084  1.00000000
Box.test.2(residuals(ajuste1),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")                   # queremos valores >0.05 en todos
##      Retard  p-value
## [1,]      6 6.49e-06
## [2,]     12 5.80e-07
## [3,]     18 2.54e-06
## [4,]     24 1.10e-07
## [5,]     30 6.80e-07
## [6,]     36 8.00e-08
## [7,]     42 4.30e-07
## [8,]     48 3.00e-07
acf(ajuste1$residuals, lag.max = 48, xlab = "Retardo", # aquí vemos que sólo se salen en multiplos de 12 y que hay patrón de simetría # para la practica es medio verdad, decidimos hacer lo otro
    main= "Función de autocorrelación simple")

pacf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

# Se propone un SARIMA(1,0,1)x(0,0,0)12 + MU
ajuste2 <- Arima(datos.train.ts,
                 order = c(1,0,1),
                 seasonal = list(order = c(0,0,0), period = 12),
                 method = "ML")

coeftest(ajuste2)
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.9947456  0.0054474 182.610 < 2.2e-16 ***
## ma1       -0.7254036  0.0434252 -16.705 < 2.2e-16 ***
## intercept  5.5238319  0.5351063  10.323 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajuste2)
##                  ar1         ma1   intercept
## ar1        1.0000000 -0.27585575 -0.01079400
## ma1       -0.2758557  1.00000000 -0.01006758
## intercept -0.0107940 -0.01006758  1.00000000
Box.test.2(residuals(ajuste2),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box") # test de jung-box
##      Retard    p-value
## [1,]      6 0.69073487
## [2,]     12 0.22973902
## [3,]     18 0.21893065
## [4,]     24 0.03960481
## [5,]     30 0.08546994
## [6,]     36 0.03580800
## [7,]     42 0.04708855
## [8,]     48 0.06789183
acf(ajuste2$residuals, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(ajuste2$residuals, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

# Se propone un SARIMA(1,0,1)x(1,0,0)12 + MU
ajuste3 <- Arima(datos.train.ts,
                 order = c(1,0,1),
                 seasonal = list(order = c(1,0,0), period = 12),
                 method = "ML")

coeftest(ajuste3)
## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## ar1        0.9922742  0.0072296 137.2515 < 2.2e-16 ***
## ma1       -0.7393339  0.0458966 -16.1087 < 2.2e-16 ***
## sar1       0.2426026  0.0734512   3.3029 0.0009569 ***
## intercept  5.5194141  0.4959124  11.1298 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Se propone un SARIMA(0,1,1)x(1,0,0)12 + MU
ajuste4 <- Arima(datos.train.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(1,0,0), period = 12),
                 method = "ML")

coeftest(ajuste4)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.743005   0.044523 -16.6881 < 2.2e-16 ***
## sar1  0.238164   0.072912   3.2664  0.001089 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajuste4)
##             ma1       sar1
## ma1   1.0000000 -0.1216682
## sar1 -0.1216682  1.0000000
Box.test.2(residuals(ajuste4),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box") # test de jung-box
##      Retard   p-value
## [1,]      6 0.7123801
## [2,]     12 0.9055981
## [3,]     18 0.8138609
## [4,]     24 0.5663195
## [5,]     30 0.7499105
## [6,]     36 0.6405039
## [7,]     42 0.6366659
## [8,]     48 0.7534805
calculoExplicativasCalendarioCaceres <- function(variableFecha, domingoYFestivosJuntos){
  
  #######################################
  #     Creación de todas las fechas    #
  #######################################
  
  if (month(max(variableFecha)) %in% c(1,3,5,7,8,10,12)) {
    diasHastaFinMes <- 30
  } else if (month(max(variableFecha)) %in% c(4,6,9,11)) {
    diasHastaFinMes <- 29
  } else if (year(max(variableFecha))%%4==0) {
    diasHastaFinMes <- 28
  } else {diasHastaFinMes <- 27}
  
  todasLasFechas <- data.frame(fechas=seq(min(variableFecha),
                                          max(variableFecha)+diasHastaFinMes,
                                          by="days"))
  
  #######################################
  #     Cálculo de la Semana Santa      #
  #######################################
  
  domingoResurrecion <- as.Date(Easter(year(min(variableFecha)):year(max(variableFecha))))
  lunesPascua <- domingoResurrecion + 1
  sabadoSanto <- domingoResurrecion - 1
  viernesSanto <- domingoResurrecion - 2
  juevesSanto <- domingoResurrecion - 3
  
  # Se unen y ordenan todos los días que forman la Semana Santa
  semanaSanta <- sort(c(juevesSanto, viernesSanto, sabadoSanto, domingoResurrecion, lunesPascua))
  
  # Se pone en formato data.frame y se añade un indicador
  semanaSanta <- data.frame(fechas=semanaSanta, semanaSanta=rep(1,length(semanaSanta)))
  
  # Se añaden a la tabla maestra de fechas
  todasLasFechas_2 <- merge(x = todasLasFechas, y = semanaSanta, by = "fechas", all.x = TRUE)
  
  # Se reemplazan los NAs por 0
  todasLasFechas_2$semanaSanta[is.na(todasLasFechas_2$semanaSanta)] <- 0
  
  
  ######################################
  #     Cálculo de la variable dt      #
  ######################################
  
  # 1. Definición de festivos:
  ############################
  
  calendario <- todasLasFechas
  
  calendario$diaSemana <- as.factor(wday(calendario$fecha))
  calendario$diaMes <- as.factor(day(calendario$fecha))
  calendario$mes <- as.factor(month(calendario$fecha))
  calendario$anyo <- as.factor(year(calendario$fecha))
  
  calendario$p_01ene <- ifelse(calendario$diaMes==1 & calendario$mes==1, 1, 0)
  calendario$p_06ene <- ifelse(calendario$diaMes==6 & calendario$mes==1, 1, 0)
  calendario$p_19mar <- ifelse(calendario$diaMes==19 & calendario$mes==3, 1, 0)
  calendario$p_01may <- ifelse(calendario$diaMes==1 & calendario$mes==5, 1, 0)
  calendario$p_15ago <- ifelse(calendario$diaMes==15 & calendario$mes==8, 1, 0)
  calendario$p_12oct <- ifelse(calendario$diaMes==12 & calendario$mes==10,1, 0)
  calendario$p_01nov <- ifelse(calendario$diaMes==1 & calendario$mes==11, 1 ,0)
  calendario$p_06dic <- ifelse(calendario$diaMes==6 & calendario$mes==12, 1 ,0)
  calendario$p_08dic <- ifelse(calendario$diaMes==8 & calendario$mes==12, 1 ,0)
  calendario$p_25dic <- ifelse(calendario$diaMes==25 & calendario$mes==12, 1 ,0)
  
  # Fiestas regionales: día de Extremadura
  calendario$p_08sep <- ifelse(calendario$diaMes==8 & calendario$mes==9, 1 ,0)
  # Fiestas locales: San Jorge
  calendario$p_23abr <- ifelse(calendario$diaMes==23 & calendario$mes==4, 1 ,0)
  
  calendario$festivo <- rowSums(subset(calendario, select=p_01ene:p_25dic))
  
  # La definición de la variable dt varía según la opción domingoYFestivosJuntos.
  
  if (domingoYFestivosJuntos==0){
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
    
    # Días laborables: todos menos sábados y domingos
    calendario$laborable <- 1-calendario$sabado-calendario$domingo
    
  } else {
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
    # Domingo=1 si domingo=1 o festivo=1
    calendario$domingo <- ifelse(calendario$domingo==1 | calendario$festivo==1, 1 ,0)
    
    # Días laborables: todos menos sábados y domingos/festivos
    calendario$laborable <- 1-calendario$sabado-calendario$domingo    
  }
  
  
  # 2. Definición de variable dt:
  ###############################
  
  # Se filtran las columnas de inter?s y se añade la Semana Santa
  
  calendario_2 <- calendario[, c("fechas", "mes", "anyo", "sabado", "domingo", "laborable", "festivo")]
  
  todasLasFechasFinal <- merge(x = todasLasFechas_2, y = calendario_2,
                               by = "fechas", all.x = TRUE)
  
  # Agregamos la serie a nivel a?o-mes
  
  calendarioAnyoMes <- aggregate(todasLasFechasFinal[,c("sabado","domingo",
                                                        "laborable", "semanaSanta", "festivo")],
                                 by=list(mes=todasLasFechasFinal$mes,
                                         anyo=todasLasFechasFinal$anyo),
                                 "sum")
  
  # Se calcula la variable dt:
  
  calendarioAnyoMes$dt <- calendarioAnyoMes$laborable-(5/2)*(calendarioAnyoMes$sabado+calendarioAnyoMes$domingo)
  
  
  ######################################
  #     Cálculo de años bisiestos      #
  ######################################
  
  calendarioAnyoMes$anyoNum <- as.numeric(levels(calendarioAnyoMes$anyo))[calendarioAnyoMes$anyo]
  
  calendarioAnyoMes$bisiesto <- ifelse(calendarioAnyoMes$mes==2 &(calendarioAnyoMes$anyoNum %% 4)==0, 1 ,0)
  
  
  #######################################################
  #     Tabla final con explicativas de calendario      #
  #######################################################
  
  if (domingoYFestivosJuntos==0){
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto", "festivo")])
  } else {
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto")])
  }
  
  return(explicativasCalendario)
  
}
explicativasCalendarioTrain <- calculoExplicativasCalendarioCaceres(datos.train$FECHA, domingoYFestivosJuntos=0)

calendarioTrain <- 
  as.matrix(
    explicativasCalendarioTrain[,c("semanaSanta", "dt", "bisiesto", "festivo")]
  )
ajuste4conCalendario <- Arima(datos.train.ts,
                              order = c(0,1,1),
                              seasonal = list(order = c(1,0,0), period = 12),
                              xreg = calendarioTrain,
                              method = "ML")

coeftest(ajuste4conCalendario)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7199518  0.0483371 -14.8944 < 2.2e-16 ***
## sar1         0.1968911  0.0733748   2.6834  0.007289 ** 
## semanaSanta -0.0115610  0.0118108  -0.9788  0.327654    
## dt           0.0150418  0.0055963   2.6878  0.007192 ** 
## bisiesto    -0.0227074  0.1049431  -0.2164  0.828693    
## festivo     -0.0614426  0.0211852  -2.9003  0.003729 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# decidimos quitar año bisiesto
calendarioTrain <- 
  as.matrix(
    explicativasCalendarioTrain[,c("semanaSanta", "dt", "festivo")]
  )

# Se podría incluir tambi?n la variable festivo

ajuste5conCalendario <- Arima(datos.train.ts,
                              order = c(0,1,1),
                              seasonal = list(order = c(1,0,0), period = 12),
                              xreg = calendarioTrain,
                              method = "ML")

coeftest(ajuste5conCalendario)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7202627  0.0482284 -14.9344 < 2.2e-16 ***
## sar1         0.1975910  0.0732627   2.6970  0.006996 ** 
## semanaSanta -0.0113950  0.0117906  -0.9664  0.333821    
## dt           0.0150369  0.0055984   2.6859  0.007233 ** 
## festivo     -0.0607887  0.0209877  -2.8964  0.003775 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajuste5conCalendario)
##                      ma1          sar1 semanaSanta            dt     festivo
## ma1          1.000000000 -0.1395694285 -0.02512203 -0.0060496156 -0.07048090
## sar1        -0.139569429  1.0000000000 -0.07245826 -0.0004302792 -0.03481278
## semanaSanta -0.025122031 -0.0724582600  1.00000000  0.0491736995  0.16027535
## dt          -0.006049616 -0.0004302792  0.04917370  1.0000000000  0.00428523
## festivo     -0.070480901 -0.0348127821  0.16027535  0.0042852305  1.00000000
Box.test.2(residuals(ajuste5conCalendario),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.7984990
## [2,]     12 0.9116931
## [3,]     18 0.8911888
## [4,]     24 0.7504270
## [5,]     30 0.9289244
## [6,]     36 0.8699718
## [7,]     42 0.8579565
## [8,]     48 0.9146346
# [7]. Identificaci?n de outliers

# Esto identifica outliers

listaOutliersTrain <- locate.outliers(ajuste5conCalendario$residuals,
                                      pars = coefs2poly(ajuste5conCalendario),
                                      types = c("AO", "LS", "TC"),
                                      cval=3)

listaOutliersTrain$abststat=abs(listaOutliersTrain$tstat)

# Cruzamos con la tabla original para obtener la fecha

datos.train$ind <- as.numeric(rownames(datos.train))
listaOutliersTrainFecha <- merge(listaOutliersTrain, datos.train[,c("ind", "FECHA")], by = "ind")

arrange(listaOutliersTrainFecha,desc(listaOutliersTrainFecha$abststat))
##   ind type    coefhat     tstat abststat      FECHA
## 1  99   LS -0.4456849 -3.135621 3.135621 2011-03-01
## 2 126   TC -0.5369467 -3.119897 3.119897 2013-06-01
## 3  88   TC  0.5255352  3.053591 3.053591 2010-04-01
outliersTrain <- outliers(c( "TC","LS","AO"), c(88, 99, 126))
outliersVariablesTrain <- outliers.effects(outliersTrain, length(ajuste5conCalendario$residuals))
calendarioMasOutliersTrain <- as.matrix(cbind(calendarioTrain, outliersVariablesTrain))
ajuste5conCalendarioYOutliers <- Arima(datos.train.ts,
                                               order = c(0,1,1),
                                               seasonal = list(order = c(1,0,0), period = 12),
                                               xreg = calendarioMasOutliersTrain,
                                               method = "ML")

coeftest(ajuste5conCalendarioYOutliers)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7705343  0.0443086 -17.3902 < 2.2e-16 ***
## sar1         0.2043685  0.0742688   2.7517  0.005928 ** 
## semanaSanta -0.0133834  0.0115593  -1.1578  0.246947    
## dt           0.0138329  0.0055408   2.4966  0.012541 *  
## festivo     -0.0611110  0.0203547  -3.0023  0.002679 ** 
## TC88         0.4155507  0.1833772   2.2661  0.023445 *  
## LS99        -0.3817126  0.1443567  -2.6442  0.008188 ** 
## AO126       -0.4052362  0.2049407  -1.9773  0.048004 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.arma(ajuste5conCalendarioYOutliers)
##                     ma1         sar1 semanaSanta           dt      festivo
## ma1          1.00000000 -0.099074246 -0.03840997  0.016422714 -0.033633401
## sar1        -0.09907425  1.000000000 -0.05525868 -0.011715920 -0.030682112
## semanaSanta -0.03840997 -0.055258676  1.00000000  0.054744004  0.152545681
## dt           0.01642271 -0.011715920  0.05474400  1.000000000  0.009458039
## festivo     -0.03363340 -0.030682112  0.15254568  0.009458039  1.000000000
## TC88         0.12700224 -0.003696872 -0.12348004 -0.011857780  0.053345083
## LS99         0.10033436  0.160271312 -0.07843167 -0.026932065  0.038140304
## AO126        0.16167463  0.024750256  0.01876841  0.137038092  0.043949322
##                     TC88        LS99      AO126
## ma1          0.127002243  0.10033436 0.16167463
## sar1        -0.003696872  0.16027131 0.02475026
## semanaSanta -0.123480045 -0.07843167 0.01876841
## dt          -0.011857780 -0.02693207 0.13703809
## festivo      0.053345083  0.03814030 0.04394932
## TC88         1.000000000  0.25658620 0.02030734
## LS99         0.256586201  1.00000000 0.02203153
## AO126        0.020307344  0.02203153 1.00000000
Box.test.2(residuals(ajuste5conCalendarioYOutliers),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.9386016
## [2,]     12 0.9353550
## [3,]     18 0.8785053
## [4,]     24 0.8525545
## [5,]     30 0.9582087
## [6,]     36 0.9668843
## [7,]     42 0.9358235
## [8,]     48 0.9096749

[7]. Predicción

explicativasCalendarioTest <- calculoExplicativasCalendarioCaceres(datos.test$FECHA,domingoYFestivosJuntos=0)
calendarioTest <- as.matrix(explicativasCalendarioTest[,c("semanaSanta", "dt", "festivo")])
# Nos valen los ?ltimos 12 datos de la variable (12 1's)
# para el valor futuro de la variable

outliersVariablesTest <- outliersVariablesTrain[181:192,]
calendarioMasOutliersTest <- as.matrix(cbind(calendarioTest, outliersVariablesTest))
calendarioMasOutliersTest
##       semanaSanta   dt festivo         TC88 LS99 AO126
##  [1,]           0  3.0       2 3.927514e-15    1     0
##  [2,]           0  0.0       0 2.749260e-15    1     0
##  [3,]           0 -4.0       1 1.924482e-15    1     0
##  [4,]           5  2.0       0 1.347137e-15    1     0
##  [5,]           0  3.0       1 9.429961e-16    1     0
##  [6,]           0 -5.0       0 6.600972e-16    1     0
##  [7,]           0  3.0       0 4.620681e-16    1     0
##  [8,]           0 -0.5       1 3.234477e-16    1     0
##  [9,]           0 -1.5       0 2.264134e-16    1     0
## [10,]           0  3.0       1 1.584893e-16    1     0
## [11,]           0 -1.5       1 1.109425e-16    1     0
## [12,]           0 -0.5       3 7.765978e-17    1     0
prediccion <- as.data.frame(predict(ajuste5conCalendarioYOutliers, 
                                    n.ahead=12,
                                    newxreg=calendarioMasOutliersTest))

#L?mites de confianza al 95%

U <- exp(prediccion$pred + 2*prediccion$se)
L <- exp(prediccion$pred - 2*prediccion$se)

datos.pred <- data.frame(FECHA = datos.test$FECHA, 
                         Prediccion = exp(prediccion$pred + 0.5 * prediccion$se^2),
                         LimSup = U, 
                         LimInf = L)
datos.real.pred <- merge(datos, datos.pred, by = "FECHA", all.x = T)

datos.real.pred$REAL <- datos.real.pred$TARGET

datos.real.pred$TARGET[datos.real.pred$FECHA>as.Date('01/12/2018',format='%d/%m/%Y')] <- NA

grafico <- ggplot(data = datos.real.pred) +
  geom_line(aes(x= FECHA, y = REAL), color = 'blue',
            alpha = 0.4, size = 0.8) +
  geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
            alpha = 0.8, size = 0.8) +
  geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
            size = 1)   +
  geom_line(aes(x= FECHA, y = LimSup), color = 'orange',
            size = 1)  +
  geom_line(aes(x= FECHA, y = LimInf), color = 'orange',
            size = 1) +
  xlab('FECHA') + ylab('Matriculaciones')

ggplotly(grafico)